molmass_glucose = 180; % Molar mass of glucose
clc;
close all;
%% Do 37C
tempmat = monod_mat_glucose_25C;
smat = size(tempmat);

conc1 = tempmat(:,1)/100*1000/molmass_glucose*1000; % Convert to mM
grmax1 = nanmean(tempmat(:,2:4),2);
grmax_err1 = nanstd(tempmat(:,2:4),0,2);

% Do fit on corrected concentration (in mM)
figure;
errorbar(conc1,grmax1,grmax_err1, 'bo', 'linewidth', 1);
    
hold on;
smat = size(tempmat);

set(gcf, 'position', [0 0 400 300])
set(gca, 'fontsize', 20);
xlabel('Concentration (mM)')
ylabel('Growth rate (1/h)')

box off;

% Monod fit
conc_tot = [conc1];
grmax_tot = [grmax1];

monodfit = @(x,s) x(1).*s./(x(2) + s) - x(3);
x0 = [1 0.01 0.01];
[x,resnorm,residual,exitflag,output,lambda,jacobian] = lsqcurvefit(monodfit, x0, conc_tot,grmax_tot);

hold on;
xfit = linspace(0,max(conc_tot), 1000);
plot(xfit, monodfit(x,xfit), 'b', 'linewidth', 1);

ci = nlparci(x,residual,'jacobian',jacobian);

err_x1 = abs(ci(1,2)-x(1));
err_x2 = abs(ci(2,2)-x(2));

disp('Fit and error of KM (mM)');
disp(x(2));
disp(err_x2);


%% Do 37C
tempmat = monod_mat_glucose_30C;
smat = size(tempmat);

conc1 = tempmat(:,1)/100*1000/molmass_glucose*1000; % Convert to mM
grmax1 = nanmean(tempmat(:,2:4),2);
grmax_err1 = nanstd(tempmat(:,2:4),0,2);

% Do fit on corrected concentration (in mM)
hold on;
errorbar(conc1,grmax1,grmax_err1, 'ko', 'linewidth', 1);
    
hold on;
smat = size(tempmat);

set(gcf, 'position', [0 0 400 300])
set(gca, 'fontsize', 20);
xlabel('Concentration (mM)')
ylabel('Growth rate (1/h)')

box off;

% Monod fit
conc_tot = [conc1];
grmax_tot = [grmax1];

monodfit = @(x,s) x(1).*s./(x(2) + s);
x0 = [1 0.01];
[x,resnorm,residual,exitflag,output,lambda,jacobian] = lsqcurvefit(monodfit, x0, conc_tot,grmax_tot);

hold on;
xfit = linspace(0,max(conc_tot), 1000);
plot(xfit, monodfit(x,xfit), 'k', 'linewidth', 1);

ci = nlparci(x,residual,'jacobian',jacobian);

err_x1 = abs(ci(1,2)-x(1));
err_x2 = abs(ci(2,2)-x(2));

disp('Fit and error of KM (mM)');
disp(x(2));
disp(err_x2);

%% Do 37C
tempmat = monod_mat_glucose_37C;
smat = size(tempmat);

conc1 = tempmat(:,1)/100*1000/molmass_glucose*1000; % Convert to mM
grmax1 = nanmean(tempmat(:,2:4),2);
grmax_err1 = nanstd(tempmat(:,2:4),0,2);

% Do fit on corrected concentration (in mM)
hold on;
errorbar(conc1,grmax1,grmax_err1, 'ro', 'linewidth', 1);
    
hold on;
smat = size(tempmat);

set(gcf, 'position', [0 0 400 300])
set(gca, 'fontsize', 20);
xlabel('Concentration (mM)')
ylabel('Growth rate (1/h)')

box off;

% Monod fit
conc_tot = [conc1];
grmax_tot = [grmax1];

monodfit = @(x,s) x(1).*s./(x(2) + s);
x0 = [1 0.01];
[x,resnorm,residual,exitflag,output,lambda,jacobian] = lsqcurvefit(monodfit, x0, conc_tot,grmax_tot);

hold on;
xfit = linspace(0,max(conc_tot), 1000);
plot(xfit, monodfit(x,xfit), 'r', 'linewidth', 1);

ci = nlparci(x,residual,'jacobian',jacobian);

err_x1 = abs(ci(1,2)-x(1));
err_x2 = abs(ci(2,2)-x(2));

disp('Fit and error of KM (mM)');
disp(x(2));
disp(err_x2);